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A time-dependent Ginzburg-Landau model of plastic deformation in two-dimensional solids is 
presented. The fundamental dynamic variables are the displacement field u and the lattice velocity 
v = du/dt. Damping is assumed to arise from the shear viscosity in the momentum equation. 
The elastic energy density is a periodic function of the shear and tetragonal strains, which enables 
formation of slips at large strains. In this work we neglect defects such as vacancies, interstitials, or 
grain boundaries. The simplest slip consists of two edge dislocations with opposite Burgers vectors. 
The formation energy of a slip is minimized if its orientation is parallel or perpendicular to the flow 
in simple shear deformation and if it makes angles of ±7r/4 with respect to the stretched direction in 
uniaxial stretching. High-density dislocations produced in plastic flow do not disappear even if the 
flow is stopped. Thus large applied strains give rise to metastable, structurally disordered states. 
We divide the elastic energy into an elastic part due to affine deformation and a defect part. The 
latter represents degree of disorder and is nearly constant in plastic flow under cyclic straining.® 

PACS numbers: 62.20.Fe, 61.72.Lk, 81.40.Lm, 83.20.Jp 



I. INTRODUCTION 

Plastic flow has long been studied in crystalline and 
amorphous solids and in glassy polymers. In crystals 
irreversible motions of dislocations give rise to plastic 
deformation and large strains produce high-density dis- 
locations [1,2]. The nonlinear flow properties are very 
sensitive to the amount of such defects and strongly de- 
pendent on the deformation history. Simulations of dy- 
namics of dislocation lines have recently been performed 
but are still most difficult [3-5]. 

In amorphous solids at low temperature T [6-12], 
salient features are asfollows. (i) Shear strains tend to 
be localized in narrow shear bands in plastic flow above 
a yield stress. The width of such shear bands is micro- 
scopic in the initial stage [13] but can grow to mesoscopic 
sizes [14], sometimes resulting in fracture. Shear bands 
were numerically realized at large shear strains in molec- 
ular dynamics (MD) simulations of two-dimensional (2D) 
two-component glasses [15,16] and in simulations of a 2D 
phenomenological stochastic model [17]. (ii) As another 
aspect, in 3D MD simulations on model two-component 
glasses at low T, Takeuchi et al. [18] observed hetero- 
geneities among mobile and immobile regions after ap- 
plication of shear strains. In 2D and 3D MD simulations 
on model supercooled binary mixtures above T g , simi- 
lar or much more extended dynamic heterogeneities have 
been detected in quiescent states [19-22] and found to 
be sensitively suppressed by applied shear flow [23] . (iii) 
Furthermore, at higher T(> 0.7 - 0.8T g ) (where T g is 
the glass transition temperature), shear deformation oc- 
curs quasi-homogeneously (still involving many particles 
in each configuration change), leading to highly viscous 
non-Newtonian behavior [23-26]. 

Glassy polymers also behave analogously [12,27,28], 
where shear bands appear above a yield stress at low 



T [29], and highly viscous non-Newtonian flow and sig- 
nificant elongation of the chain shapes occur at elevated 
T [24,30]. In glassy polymers, the entropic stress arising 
from molecular orientations becomes significant at large 
strains and such systems behave like cross-linked rubbers 
[12,27,28]. 

Recently much attention has been paid to jamming 
rheology observed in sheared states of supercooled liq- 
uids, soft glassy materials such as dense microemulsions, 
or granular materials [31]. In these systems, the ther- 
mal agitation effect is very small if the particle size is 
large, but universal constrained dynamics is realized un- 
der external forces. In supercooled liquids (at relatively 
high T) [23] and dense microemulsions (at effectively low 
T) [32,33], mesoscopic dynamic heterogeneity and strong 
shear-thinning behavior have been observed, but shear 
bands have not been identified. In granular materials (at 
effectively zero temperature), strain localization is most 
conspicuous [34-36]. 

In our recent work [37] we have constructed 2D nonlin- 
ear strain theory taking into account the underlying local 
periodic lattice structure, where the elastic energy den- 
sity is periodic with respect to the shear and tetragonal 
strains. Then in our model plastic flow starts with ap- 
pearance of slips. The simplest slips consist of two edge 
dislocations having opposite Burgers vectors with size a 
(a being the lattice constant) and they grow into meso- 
scopic shear bands as the applied strain is increased. Un- 
der uniaxial stretching [38], well developed shear bands 
were already numerically realized [16,17] and will be re- 
alized also in our simulations. These shear bands make 
angles of ±7r/4 with respect to the stretched direction in 
agreement with observations in various amorphous mate- 
rials [7-10,27] including granular materials [34,35]. These 
angles will be shown to minimize the elastic energy of in- 
cipient slips in this paper. 
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In crystalline solids, dislocation pairs in 2D or disloca- 
tion loops in 3D forming slip lines or surfaces should be 
nucleated at the inception of plastic deformation (in ad- 
dition to preexisting dislocations). In amorphous solids, 
it has been controversial on whether dislocations them- 
selves can be well-defined or not [6,9]. It is not obvious 
how to characterize the local rearrangement processes as 
to their shapes and sizes. In MD simulations on binary 
mixtures, they have been detected as cluster-like objects 
with various visualization methods [18-24,39]. If the 
size-ratio of the constituent two species is chosen such 
that crystallization is most suppressed, strong frustra- 
tion occurs in the packing of large and small particles 
in jammed states. Then the local crystalline order can 
be well-defined only over short distances. However, in a 
2D amorphous soap bubble raft, Argon and Kuo [40] ob- 
served that nuclcation of a dislocation pair gave rise to a 
small-scale slip but such dislocations did not glide more 
than 2-3 bubble distances. It is worth noting that Deng 
et al. [15] found extended slip-like strain localization in 
2D MD simulations with the size-ratio rather close to 1 
(see comment (iv) in the last section). Notice that the 
displacement field around a slip is localized because the 
two constituent dislocations have opposite Burgers vec- 
tors. As a result, small-scale slips should be well-defined 
even in amorphous solids as long as the slip size does not 
much exceed the range of the local crystal structure. 

The plastic flow phenomena are thus very complex, be- 
ing influenced by many factors, but they are ubiquitous 
in various kinds of solid-like materials. The purpose of 
this paper is to present a well-defined Ginzburg-Landau 
model consistently taking account of nonlinear elasticity. 
A merit of this approach is that we can put emphasis on 
any aspects of the phenomena by controlling the param- 
eters or changing the model itself. We will examine (i) 
the fundamental flow units, slips, in detail numerically 
and analytically and (ii) plastic flow numerically in sim- 
ple shear and clongational (stretching) deformation. To 
make this paper simplest, as it is the first detailed ex- 
position of our scheme, we will neglect (i) vacancies and 
interstitials, or a variable m representing the local free- 
volume. A dynamic model including such an additional 
degree of freedom has already been presented in our pre- 
vious work [37]. We will also neglect (ii) the configu- 
rational frustration effect induced by the size difference 
between the two species. Introducing these two ingredi- 
ents will constiture future development of our scheme. 

This paper is organized as follows. In Section 2 we will 
present our dynamic model and explain our numerical 
method. In Section 3 we will discuss the simplest form 
of slips numerically obtained from our nonlinear elastic- 
ity theory. We will also derive some analytic expressions 
for the slip formation energy under general strain field 
starting with the Peach-Koehler theory [41] and compare 
them with numerical results. In our scheme stationary 
slip solutions exist for small externally applied strains, 



where the force balance is achieved in the presence of 
the Peierls potential energy for the dislocation position 
in crystals [42,43]. Section 4 will present numerical re- 
sults of the stress-strain relations and the patterns of the 
strains, the elastic energy density, and the displacement 
vector under applied strains. We will also give a method 
of dividing the elastic energy into an elastic part due to 
affinc deformation and a defect part. 



II. MODEL EQUATIONS 

A. Elastic energy 

Recently Doi et al. [44] described plastic flow in a 
highly viscous 2D crystal phase of block copolymers in 
the framework of a ID Frenkcl-Kontrowa model [45] as- 
suming ID sliding motions. Hereafter we present a 2D 
Frenkel-Kontrowa model to describe plastic flow in 2D 
[37]. In terms of the displacement vector u = (u x ,u y ) 
from a reference crystal state, we define the strain com- 
ponents as 

ei = V x u x + \7 y Uy, 
e 2 = V x u x - VyUy, 

e?, = V x Uy + VyU x , (2.1) 

where \7 X = d/dx and \7 y = d/dy. We call e\ the di- 
lation strain, e 2 the tetragonal strain, and e 3 the shear 
strain. If we suppose a 2D triangular lattice with lattice 
constant a, the elastic energy should be invariant with 
respect to the rotations of the reference frame by ±mr/3 
(n = 1, 2, • • •). Due to this symmetry, the elasticity must 
be isotropic in the harmonic approximation [46], being 
characterized by the bulk and shear moduli, K and /xo, 
but it depends on the oricntational angle 9 of one of the 
crystal axes with respect to the x axis for large shear 
strains. Under rotation of the reference frame by 9, the 
shear strains e 2 and e 3 are changed to e' 2 and e 3 , where 
[46] 

e' 2 = e 2 cos 29 + e 3 sin 29, 

e' 3 = e 3 cos 29 - e 2 sin 29. (2.2) 

These relations are obtained from the orthogonal trans- 
formations, r' =JJ -r and u' =[] u, where r' = (x',y') 
and u 1 = (u' x , u' y ) represent the position and the displacc- 

ment, respectively, in the new reference frame, and U— 
{Uij} with U xx = U yy = cos9 and U xy = —U yx = sin#. 
The e' 2 and e 3 are the tetragonal and shear strains, re- 
spectively, in the new reference frame where the x' axis 
is along one of the crystal axes. 

The elastic energy is written as F e i = J drf c \ with the 
elastic energy density in the form, 



which is independent of the rotation strain, 

W = V X Uy - VyU X . 



(2.3) 



(2.4) 



Note that ei and lo are invariant with respect to the ro- 
tation of the reference frame. The simplest form of $ is 
given by 



3 — cos ir(V3e 3 — e' 2 ) 



— cos7r(\/3e 3 + e' 2 ) — cos(27re 2 ) 



(2.5) 



This function is invariant with respect to the rotation 
9 — > 9 + 7r/3, is a periodic function of e' 3 with period 
2 / \/3 for e' 2 = (simple shear deformation) , and becomes 
(e| + e§)/2 for small strains. In Fig.l we display §(e 3 , e 2 ) 
supposing 9 = 0, where one of the crystal axes is along 
the x axis. We can see a hexagonal lattice structure in the 
63-62 plane. As a characteristic feature, Fig. 2 shows that 
it is almost isotropic or is a function of e = (e§ + e^) 1 ^ 2 
only even for rather large e (< 0.5) (in the unit cell at 
the origin). In fact, we have the Taylor expansion, 



$(e3,e 2 ) 



-e — 



-e 1 + —e° 
+ 72 



+ ^(et-et)(e i -16ete 2 3 ) + 0(e*) 



(2.6) 



If we set e 3 = ecosx and e 2 = esinx, the fourth term is 
rewritten as — (7r 4 /720)e 6 cos(6x) and is known to keep 
the invariance x ~ > X + it/3. Anisotropy appears in the 
terms of order e 6 , but the anisotropic fourth term is at 
most 10% of the isotropic third term. 

We examine elastic stability of homogeneously strained 
states. That is, we superimpose infinitesimal strains, <5e3 
and Se 2 , on e 3 and e 2 assumed to be homogeneous. The 
second order terms in $(e3 + 5e 3 , e 2 + Se 2 ) read 

5(2)$ = -^(Se 3 ) 2 + $ 23 Se 3 5e 2 + ^f{5e 2 ) 2 , (2.7) 

where $ a/ 3 = d 2 &/de a dep (a,/3 = 2,3). In the stable 
regions the above second order contribution should be 
positive-definite or the two eigenvalues, Ai and X 2 , of 
the 2x2 matrix {Qap} should be both positive. In the 
appendix we will derive this linear stability criterion by 
solving the linearized version of our dynamic model to 
be presented in the next section. Fig. 3 shows that the 
system is stable for e < 0.3, where $ depends almost 
only on e. This result readily follows if we neglect the 
anisotropy in $ by setting <J>(e 2 , e 3 ) = G(e 2 ). Then some 
calculations yield 



$ a/3 = 2G' ( 5 Q(3 + 4G"e Q e / 5, 



(2.8) 



where G' = dG(e 2 )/de 2 and G" = dG'(e 2 )/de 2 . From 
(2.6) we find G' ^ 1/2 and G" = -ir 2 /4. The deter- 
minant of {$ Q/3 } becomes 4G'(G' + 2G"e 2 ). Thus the 
stability condition becomes 



e < (G7|2G"|) 1/2 ^0.3, 



(2.9) 



around the origin in the e 3 -e 2 plane. The elastic instabil- 
ity with a negative eigenvalue (Ai < or \ 2 < 0) causes 
rapid relaxation processes resulting in localized slips and 
stress release. In plastic flow the stability condition is 
mostly satisfied throughout the system (see Fig. 11). 

If we assume that small solid elements are rotated 
without shape changes with the local angular velocity 
(dvy/dx — dv x /dy)/2, the orientation angle 9 is related 
to the rotation strain as 



9 = -w + 6 , 



(2.10) 



where 9 is the initial value independent of time t. This 
relation was assumed in our previous paper [37]. How- 
ever, it becomes not well defined when the three crystal 
axes are rotated differently For example, let a crystal 
with the three axes at 9 = 0,7r/3, and — n/3 be affmcly 
deformed by a simple shear deformation given by u x = 
and u y = (e 3 = 7 and e 2 = 0). Then the the first axis 
along the x direction is not rotated, while the other axes 
are rotated by 



se± = ± 



tan 



V3 



1± V3 7 



(2.11) 



For 7 < 1, we have S9± = -37/4, so (2.10) holds for the 
average of the three rotation angles, (0 + S9 + + S9- )/3 = 
—7/2 = lu/2. Thus there remains ambiguity in (2.10) 
particularly at large strains e > 0.5 [47]. In this paper, 
in view of the virtual isotropic behavior of our elastic 
energy for e < 0.5, we will simply set 



= 0, 



(2.12) 



throughout the system. To check the appropriateness 
of this assumption, we have also performed simulations 
with various homogeneous 9 held fixed. The resultant 
patterns and stress-strain relations are insensitive to 9 
and very similar to those obtained for 9 = to follow. 



B. Dynamic equations 

The elastic stress tensor should be defined for general 
strains. Let us change an arbitrary displacement u by an 
infinitesimal amount Su as u — > u + Su. The incremental 
change of the elastic energy density is written as [46] 



d 



(2.13) 



This is just the definition of the elastic stress tensor cry. 
Under (2.10) [47] it is symmetric and its components are 
expressed as 



&xx = K ei + poA 2 , Uyy 



K ei - p A 2 , 



(2.14) 



where A a = <9$(e 3 , e 2 )/de a are written as 



3tt 



coslV^nes) sin(7re2) + sin(27re2) 



A 3 = j_ sin(\/37re 3 ) cos(-7re2). 

V37T 

For small strains we have A 2 = e 2 and A 3 
duce the isotropic, linear elastic theory. 
We assume that the lattice velocity 



(2.15) 



d_ 

dt 



u 



e 3 to rcpro- 



(2.16) 



obeys 



p—v = V- a + + ?7 V v + V- <r R , 



(2.17) 



where we introduce the shear viscosity 770 but neglect the 
bulk viscosity [48,49]. The <7 R = Wfj} is a symmetric 
random stress tensor and satisfies a^ x +cr^ y = 0, because 
the bulk viscosity is neglected, and is related to 770 by 
[46] 

(a? j (r,t)a? j (r',t')) = 2k B T m S(r-r')5(t-t'). (2.18) 

The mass density p will be treated as a constant in (2.17). 
This is justified when the deviation dp = p — (p) is as- 
sumed to be much smaller than the average (p) [50]. In 
the usual linear elasticity theory [48] it is related to the 
dilation strain e\ as dp = —pe\, so we are assuming 
|ei| <C 1 (and coincidence of the lattice and mass ve- 
locities [37]). In our simulation in the plastic flow regime 
at 7 — 10~ 3 , for example, |ei| attains a maximum in a 
range of 0.2-0.3 and the variance y 7 (ef) increases up to 
about 0.06. 

Due to the presence of the random stress, u and v arc 
random variables and (2.16) and (2.17) constitute non- 
linear Langevin equations [46]. Their equilibrium distri- 
bution attained in the unstrained condition is given by 
Poc exp(— F/ksT), where the total free energy F is the 
sum of the elastic energy F e i and the kinetic energy as 



"/*( 



r 1 P 2 
/el + ~ 2 V 



(2.19) 



Furthermore, if the random stress is omitted in (2.17) 
and the dynamic equations are treated as deterministic 
ones, the time derivative of F is nonnegative-definite in 
the unstrained condition as 



J t F = ~ J drJ2vo(V l v J ) 2 <0, (2.20) 



where use is made of the relation, 



ou 



(2.21) 



The stationary condition dF/dt = is attained for v = 
and V- 0= 0. This is a condition to guarantee self- 
consistency of the dynamic equations which have stable 
equilibrium solutions. 

In the appendix we will examine the linearized dy- 
namic equations of our model around homogeneously 
strained states to obtain two sound modes in the stable 
region. 



C. Dimensionless forms 

We make our equations dimensionless measuring space 
in units of the lattice constant a and time in units of 



to 



(P/Mo) 



1/2 



a. 



(2.22) 



The stress components and the elastic energy density are 
measured in units of po and the clastic energy in units of 
Pod 2 , while the strains remain unsealed. To avoid to in- 
troduce too many symbols, we will rewrite the scaled po- 
sition vector a _1 r, time r ' 1 t, displacement vector aT 1 u, 
and velocity Toa~ 1 v simply as r, t, u, and v using the 
original notation. Then, in the dimensionless a, p is 
replaced by 1 and Kq by the ratio, 



A = Ko/po- 



(2.23) 



In terms of the scaled quantities the equilibrium distri- 
bution is written as 



P (x exp 



eth J 



1 



drl -ef + $(e 3 ,e 2 ) + -v 



where 



e t h = k B T/p a 2 



(2.24) 



(2.25) 



The parameter eth represents the degree of the thermal 
fluctuations (being proportional to T) and is an impor- 
tant parameter, for example, in describing the decay of 
metastable states by thermal agitations. If the dynamic 
equation (2.17) is made dimensionless, the dimensionless 
viscosity is given by 

V * = Vo a- 1 (pp Q )- 1 / 2 . (2.26) 

In (2.18) the noise strength fce^ijo is replaced by ethVo- 



D. Numerical method 

We integrate (2.16) and (2.17) in the dimcnsionless 
units on a 128 x 128 square lattice represented by (n, m) 
(1 < n, m < 128) with x — nAx and y = mAx. The 
modulus ratio A in (2.23) is set equal to 1. For simplic- 
ity, the mesh size Ax and the dimcnsionless viscosity in 
(2.26) are set equal to 1: 

Ax = l, % * = 1. (2.27) 

The hrst relation means that the mesh size is just equal 
to the lattice constant, and the second one is rewritten 
as rj = a(pno) 1 / 2 . 

We give our boundary conditions employed. At the 
bottom y — 0, we set u = 0. At the top y = Lq = 128, 
we set u x = "/Lq and u y — in the presence of ap- 
plied shear strain 7, and we set u x = —u y = cLq/2 
in the presence of applied tetragonal strain e. In the 
x direction, we impose the periodic boundary condition, 
u(x + LQ,y,t) = u(x,y,t). The unstrained condition is 
realized for 7 = and e = 0. 

We are interested in slips across which the atomic dis- 
placement is discontinuous by the lattice constant a. In 
this paper, by setting Ax = a, we try to realize such 
singular objects numerically. For this purpose it is con- 
venient to define the strains and tensors on the middle 
points (71 + 1/2, 777 + 1/2), while the vectors® are defined 
at the lattice points (77,777). For a vector component A 
(say, A = u x ), \7 X A and S7 y A at (n+1/2, m+1/2) are de- 
fined as [A(n+1, 771+ 1) — A(n, m)+A(n+l, m)—A(n, m+ 
l)]/2 and [A(n + 1, m + 1) - A(n, m) - A(n + 1, m) + 
A(n,m + l)]/2, respectively, using A at the four points 
(77 + 1/2 ± 1/2, m + 1/2 ± 1/2). In the same manner, we 
may construct the vector V- <J at (77, m) using the stress 
components at the four points (77+ 1/2,777+ 1/2). With 
this space discretization, slips consisting of a straight line 
segment become well-defined if their angle with respect 
to the x or y axis is or 7r/4. For other slip orientations, 
zigzag points appear along the slip line segment and an 
extra elastic energy becomes needed. On the other hand, 
if we define all the quantities on (77,777), slip discontinu- 
ity takes place over a few lattice sizes, but macroscopic 
features such as the stress-strain relation remain almost 
unchanged. Furthermore, we will suppose simple shear or 
uniaxial deformation and, as will be shown in the next 
section, the preferred slip orientation angle is or 7r/4 
with respect to the x or y axis. By this reason our sim- 
ple numerical scheme seems to be allowable at least in 
this first attempt. 



III. SLIPS 



They are analogous to quantum vortex rings in superfluid 
helium [46]. Their elastic structure far from the disloca- 
tion cores may well be described by the linear elasticity 
theory, but nonlinear elasticity theory is needed (i) to 
suppress the divergence of the stress at the cores and (ii) 
to stabilize the slips themselves when they adjust to the 
crystal structure. Slips are not in a stationary state in 
the linear elasticity theory in the absence of impurities 
etc. which can trap dislocations, as will be evident in 
(3.20). In our nonlinear theory, those along the x axis 
(9 = 0) can be in a stationary metastable state if their 
length is a multiple of the lattice constant. This is consis- 
tent with the Peierls-Nabarro theory [42,43], which takes 
into account the discreteness of the crystal structure and 
gives a periodic Peierls potential energy for the position 
of the dislocation center. 



A. Slips in linear elasticity theory 

To begin with, let us write the solution of an edge dislo- 
cation as u = bu Lc = {bu)f, buy ), whose Burgers vector 
is assumed to be along the x direction and is written as 
b = (b,0). The linear elasticity theory [48] gives 



2tt 



tan 

-1 



1 



xy 



2(1 - v) x 2 +y 2 _ 



4tt(1 - v) 



(1 - 2v) In yjx 2 +y 2 + 



x 2 + y 2 



(3.1) 



In the 3D theory the dislocation line is along the z axis 
and v is 3D Poisson's ratio. In our 2D theory v is related 
to A in (2.23) as 



1 

V= 2 



1 

2A' 



(3.2) 



In our simulations A = 1 and hence v = 0. In the linear 
theory a slip is a superposition of two edge dislocations 
with opposite Burgers vectors expressed as 



M LS± = ± 



u Lc (x- -, y )-u^(x+-,y) 



(3.3) 



where i is the slip length and b = ±1 (both in units of a). 
The slip line segment is between the two points {—1/2, 0) 
and (€/2, 0). The + sign corresponds to a clockwise slip 
(type C), and the — sign to a counterclockwise slip (type 
CC). The displacement vector around a slip is clockwise 
(counterclockwise) for type C (type CC) (as will be evi- 
dent in Fig. 5 below). Across the line segment |x| < IJ2 
and y = 0, the displacement is discontinuous as 



In our model system fundamental flow units in plastic 
deformation are slips composed of a pair of edge disloca- 
tions with opposite Burger vectors (dislocation dipoles). 



,Ls±/ 



(x,y + 0)-u L x s± (x,y-0) = ±l. 
The corresponding strains are written as 



(3.4) 
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,Ls± 



„Ls± 



±{2v- 1) 
2tt(1 - v) 

tt(1 - i/) 

±1 
2tt(1 - iTj 



(x2 + y 2 ) 2 (x 2 _ + y 2 ) 2 
x + (x\ — y 2 ) X-{x 2 _— y 2 ) 

2^2 (a;2_ + y 2)2 



(4 + y 2 ) 
±<%)e(^ 2 /4-z 2 ), 



(3.5) 



where x± = x £/2, and 0(C) is the step function being 
equal to 1 for ( > and to for £ < 0. There appears 
no 5-function in the dilation and tetragonal strains. The 
strains diverge at the cores where the linear elasticity 
theory breaks down. 



B. Slips in nonlinear elasticity theory 

Next we numerically construct the corresponding slip 
solution in our nonlinear elasticity theory for integer £. 
That is, by starting with the linear solution m Ls± in (3.3) 
at t = and neglecting the random stress a r, we inte- 
grate (2.16) and (2.17) to seek the steady solution u s± 
attained after transient relaxation. That is, 



u (x, y) = lim u(x,y,t), 



(3.6) 



where u(x,y,0) = u Ls± (x,y). The core points of the 
linear solution, where the strains diverge, are placed at 
middle points (n+1/2, m+1/2) of the mesh of integration 
at t = 0. The limit u s± is a steady metastable solution 
satisfying the mechanical equilibrium condition, 



V- (7= 0. 



(3.7) 



It nearly coincides with the linear solution it Ls± far from 
the dislocation cores (\x±£/2\ 2 + y 2 > 3 in our case) and 
keeps to satisfy the slip condition (3.4). The stain and 
stress components calculated from it s± are finite even in 
the core regions. In fact, in the presence of a slip along 
the x direction, the maximum values attained by |ei|, 
|e 2 |, |e 3 |, and \<r xy \ in the core regions are about 0.18, 
0.08, 1.1, and 0.1, respectively. In Fig. 4 we show the 
x component u x = of a clockwise slip with I = 20 
in the unstrained condition, which is discontinuous by 1 
across the slip segment. 

The elastic energy of a slip in our nonlinear theory is 
then of central importance. We will neglect the Peierls 
potential energy for the time being. We generally assume 
that a slip line is oblique to the x axis making an angle 
of <p and is under externally applied strains, 



(e3> = 1, (e 2 ) 



(3.8) 



space average. The slip energy to create a single slip is 
defined by 



"slip 



F-F 



J dr(f cl 



/S). 



(3.9) 



where / e i is the elastic energy density in our nonlin- 
ear theory with one slip calculated numerically, and 
/el = ^(T' 6 ) 1S that in a homogeneously strained state. 
In the unstrained condition (7 = 6 = 0), the expression 
-Fsiip = ln£/27r(l — v) (in units of fioa 2 ) is well known 
[51]. 

Dislocation motions perpendicular to the slip line 
(climb motions) create a large number of defects (oc £) 
and the energy needed is very large, so they may be 
neglected [52]. On the other hand, dislocation motions 
along the slip line (glide motions) involve displacements 
of a relatively small number of particles (of order 10 for 
a unit-length motion as can be seen in the inset of Fig. 9) 
and hence play a major role in plastic flow. A force 
T = (T x ,T y ) acting on a dislocation under an applied 
stress {c° x } may be calculated using the Peach-Koehler 
theory [1,2,41]. For b = (b, 0) along the x axis, the com- 
ponents of the force are written as 



(3.10) 



(i) For a slip along the x axis, T y is canceled by the force 
due to an additional elastic deformation in the surround- 
ing medium. If the dislocation on the left is fixed and 
that on the right is moved by 81 along the x axis, the 
change of F s ii p is equal to T x 8l. Here the Peierls force is 
neglected. Therefore, 



d 

—F cl = T X = -oT y b. 



(3.11) 



The a X y consists of the stress produced by the dislocation 
on the left and the externally applied stress. For small 7 
we obtain 



slip 



2tt(1 - v) 



Til, 



(3.12) 



where — is for type C and + is for type CO (ii) For 
general angle ip of the slip with respect to the x axis, we 
rotate the reference frame by ip to obtain the shear strain 
7' = 7cos2(/9 — es'm2ip in the new reference frame from 
(2.2). Therefore, 



F f 



ln£ 



slip 



2tt(1 - v) 



=F (7 cos 2p — e sin 2ip)£. (3.13) 



If ip is varied in (3.13), F s i; p is minimized for 



tp=-(nn-a) (n = 0, 1, 2, • • •), 



(3.14) 



We assume (ei) = 0. The average displacement is given 
by (u x ) = fy + ex/2 and (u y ) = —ey/2. Here (• • •) is the 



where n is even for type C and odd for type CO, and a 
is determined by 



cos a = 



sin a 



(3.15) (dSui/dx-j 



The deviation Siij consists of 



For simple shear deformation with 7 > and e = 0, the 
most favorable slip orientation with the lowest F s ii p is 
<p = for type C and <p = n/2 for type CC. For uniaxial 
stretching with 7 = and e > 0, it is given by <p — — 7r/4 
for type C and (p — ir/4 for type CC, in agreement with 
the experiments [7-10,27,34,35]. 

As will be shown in the appendix, the slip directions 
determined by (3.14) are perpendicular to the wave vec- 
tors which minimize the angle-dependent sound velocity 
following from our linearized dynamic equations, or per- 
pendicular to the softest directions for the sound modes. 
This coincidence is obtained with the aid of the isotropic 
behavior (2.8), so it is not a general result for anisotropic 
solids. We remark that previous theories of strain local- 
ization [35,53,54] are based on ID analysis, where all the 
quantities vary only in one direction n normal to the 
plane of the band, and reduce to linear stability analysis 
for small amplitude perturbations for the wave vector in 
the direction of n. 

In Fig. 5 we display the displacement vector u around 
type C and type CC slips calculated in the unstrained 
condition. The slips are oriented in the most favorable 
directions in shear deformation in (a) and in uniaxial 
stretching in (b). Away from the slips the directions 
of u continuously change to those of the macroscopic 
deformation supposed to be applied. In Fig. 6 the slip 
energy F s u p of type C slips with £ = 10, 20, and 30 is 
shown as a function of the applied shear strain 7. For 
I7I < 0.05 we confirm (3.12) with the — sign. For larger 
7 the linear relation dF s u p /dj =const. docs not hold. 
For 7 < 7ci(~ —0.1) or for 7 > 7 C 2(^ 0.1), a steady 
metastable solution becomes nonexistent and, as a re- 
sult, the slip grows up to the system length or shrinks 
to vanish in the simulation. In Fig. 6 the slip energy 
/el — /ei with /g] = $(7,0) is displayed for the three 
values 7 = 0,0.065, and —0.04. Interestingly, the elas- 
tic energy density in the middle region between the two 
dislocations at the ends is decreased for positive 7 and 
increased for negative 7, giving rise to the contribution 
— £-f in the slip energy. In the core region, f e \ has two 
peaks and is rather insensitive to applied strains in our 
model. 

Furthermore, in Fig. 8 we numerically demonstrate co- 
incidence of dF s ii p /d"f and the space integral of o~ xy —j for 
a single isolated slip in simple shear deformation. This 
relation may be obtained from (2.1). If 7 is increased by 
an infinitesimal amount £7, the change of F s ii p is written 
in terms of the incremental displacement dm as 



16) 



where {of,-} is the stress in the homogeneous state 
and use is made of the relation for the space average 



the applied displacement change SixySj and the induced 
deviation <5u- localized near the slip. However, the con- 
tribution from Su'i vanishes in (3.16) from the mechanical 
equilibrium condition (3.7). Thus, 

JUii P (^7) = J dr[v xy -*%,]. (3.17) 

In Fig. 8 7 is in the range [7! < 0.1, so a xy = 7 holds 
excellently for our clastic energy Similarly, the counter- 
part of (3.17) in the case of uniaxial stretching is written 
as 

J^iip^, e) - \ J dr[a xx - a yy - a° xx + a° yy }. (3.18) 

The relations (3.17) and (3.18) hold for any strain am- 
plitudes as long as a steady slip solution can be obtained 
from (3.6), while (3.12) or (3.13) is valid only for very 
small strains (< 0.05 in our case). 



C. Peierls potential energy 

We continue to consider a slip along the x axis under 
simple shear deformation, but the slip length £ here can 
be noninteger. For general £ we modify (3.12) as 



In ^ 



slip 



2tt(1 - 



T7* + fpN(*), 



(3.19) 



where Upn(£) represents the Peierls potential energy be- 
ing zero for integer £ and positive for noninteger I. The 
force acting on the slip along the glide direction is then 
given by 

^F slip = [2tt(1 - V )\- x \ T 7 + J^pnW- (3.20) 

This force vanishes for a stationary slip. In the realm 
of linear elasticity theory, where the crystal structure is 
smoothed out, the Peierls potential energy is nonexistent 
and we are led to the following conclusion: A type CC 
slip shrinks for any £(> 1), but a type C slip shrinks for 
£ < l\ and expands for I > £\. Here we assume 7 > 
and the critical length £^ in the linear elasticity theory 
is determined by 



£\ = 1/[2tt(1 - v)i\. 



(3.21) 



If £ is fixed, we obtain a critical shear strain 7^ = 
l/[27r(l — v)£\. In our nonlinear theory, however, this 
critical length (or strain) loses its physical relevance. 

In our simulations (without the random stress) slips 
can be stationary suggesting the existence of the Peierls 
potential energy. To show this, we numerically create 
two type C slips along the x axis obtained in the limit 



(3.6); one is in the range < x < 20 with length 20, 
and the other is in the range — 1 < x < 20 with length 
21. Here the positions of the dislocation core on the left 
are different by 1 but those on the right coincide. As 
a result, the corresponding displacements M20 and U21 
are different only near the core region on the left, as can 
be seen in the inset of Fig. 9. Note that the difference 
u' = U21 — U20 is the displacement realized when the 
shorter slip grows into the longer one. Now, we calculate 
F for the interpolated displacement, 

ui = (1 - a)u 2 a + au 2 \ = u 2 o + au'. (3.22) 

As the slip energy at £ = 20 + a (0 < a < 1), 
F sVlp (£) = F - F is determined as in (3.9). The F sVlp (£) 
here depends on the displacement path connecting u 2 o 
and U21. 

In Fig. 9 the resultant energy difference Ai* 1 ,,^ = 
F s ii P (£) — -Psii P (20) is shown. From this calculation the 
Peierls potential cannot be determined uniquely but, if 
we define 

U PN (i) = F slip (£) - (1 - a)F slip (20) - aF slip (21), (3.23) 

we can know its behavior in the range < a = £ — 20 < 1 
directly from Fig. 9. We recognize that F s u p (£) takes a 
maximum at £ — ^ max between 20 and 21 and takes lo- 
cal minima at £ — 20 and 21. It follows that the stable 
force-balance condition dF s u p /£ — with d 2 F s u p /£ 2 > 
holds at £ = 20 and 21. The maximum of C/pn is about 
0.03 for 7 = 0. Fig. 9 also indicates that £ max — > 21 as 
7 -> 7d(= -0.09) and £ max -» 20 as 7 -» 7c2 (^ 0.12). 
For integer I and noninteger £', the critical strains, 7ci(^) 
and "fc2(£), with respect to shrinkage and expansion sat- 
isfy the following: (i) dF shp (£')/d£' > for any £' smaller 
than £ for 7 < 7^), (ii) dF slip (£')/d£' < for any £' 
larger than £ for 7 > "f c2 (£), and (iii) F s u p (£) is locally 
minimum at integer length £ for J c i{£) < 7 < lc2(£)- 
These are consistent with the positions of the instability 
points in Fig. 6. 



IV. PLASTIC FLOW 



In this section we induce deformation at a constant 
strain rate, 7 or e, or cyclic shear deformation for t > 
in the presence of the random stress tensor ctr with 
eth = 0.1 (except for the curve (b) in Fig. 22 where 
e t h = 0.25). At t = 0, the values of v at the lattice points 
are Gaussian random numbers with variance 0.01. The 
shear stress and the normal stress difference in the follow- 
ing figures are the space averages, (<r xy ) and (<r xx — <J yy ), 
respectively. 



A. Shear deformation 

In Fig. 10 we show the stress-strain curves obtained by 
integration of (2.16) and (2.17) under a constant shear 
rate 7. At small -f(t) = jt these curves first follow the 
curve in the homogeneous case, 

<?x V = 4$ (7) = ~y=^; sin(V37r7), (4.1) 

where A3 was introduced in (2.15). For the curve of 
7 = 10~ 4 and that of 7 = 10~ 3 (with the higher peak) 
we set u — at t — (supposing a perfect crystal) . They 
approach the elastic instability point 7 = y/3/6 = 0.289, 
where da xy /d^f tends to vanish and softening with re- 
spect to further shear deformation occurs. Then the 
shear stress drops sharply after the peak with catas- 
trophic formation of slips. Sec the two snapshots of 
^ e 3 = e3 — 7 in the figure. The orientations of the slips 
are the most favorable ones with (p — and 7r/2 as de- 
termined in (3.14). For 7 = 10~ 4 the slip formation is 
triggered at a smaller strain (= 0.25) and the spacing 
between the slips is a few times wider than in the case of 
7 = 10~ 3 . In our simulation the slip spacing depends on 
the shear rate 7 in the plastic flow regime with the other 
dimensionless parameters held fixed. For the other curve 
of 7 = 10~ 3 (with the lower peak) we put four slips with 
length 20 at t — 0, as will be illustrated in Fig. 13 below 
in detail. This curve indicates that the overshoot in the 
stress-strain relation is weakened by the initially preex- 
isting defects. In Fig. 11 we display mechanically unstable 
points (dots) for 7 = 0.29, 0.30, and 0.40 at 7 = 10~ 3 , 
where at least one of the eigenvalues Ai and A 2 of the 
matrix {$0,3} defined below (2.7) is negative. We rec- 
ognize that the stability condition (Ai > and A 2 > 0) 
are satisfied at most points in plastic flow. Fig. 12 dis- 
plays a snapshot of the displacement vector u in a 1/4 
region of the total system at 7 = 0.4 with 7 = 10~ 4 . We 
can see a number of slips and bands (aggregates of slips 
here), where ip = for type C and ip — n/2 for type 
CC. The large horizontal shear band in the lower part is 
particularly conspicuous, where the band thickness and 
the discontinuity of u x across it are both increased up to 
about 3-4. We recognize that elementary slips (dipolcs 
of edge dislocations) tend to be created around preex- 
isting ones, yielding thicker shear bands. Similar thick 
shear bands have been observed in previous simulations 
[16,17]. 

In Fig. 13 we follow time-development of 8ez and <5/ c i = 
/el — (fd) at 7 = 10~ 3 in the presence of four slips at 
t = 0, where the favorable slips (in black in the upper fig- 
ure at 7 = 0) grow and the unfavorable ones (in white) 
shrink as 7 increases. In the lower snapshots the elas- 
tic energy density deviation Sf e \ = f e \ — (f e \) is shown, 
where the black dots represent the dislocation cores. At 
7 = 0.176 we can see that the energy density between the 



two dislocations at the slip ends is decreased for the fa- 
vorable slips and is increased for the unfavorable slips, in 
accord with (3.13) and Fig. 7. At the plastic flow regime 
7 = 0.387 the initial favorable slips grow into thick layers 
where the dislocation density is very high. 

Next we apply a cyclic shear deformation, where j(t) = 
10~ 3 in the time regions nt p < t < (n+l/2)i p and j(t) = 
— 10~ 3 in the time regions (n + l/2)t p < t < (n + l)t p . 
We choose t p = 1000, so the maximum of ^(t) is 0.5 and 
the minimum is 0. For the first two cycles (n = and 1), 
Fig. 15 shows the stress-strain curve, while Fig. 16 shows 
the average elastic energy density (/ c i) with a snapshot 
of / e i at the point A where the average stress vanishes. 
As salient features, we notice (i) residual strain at van- 
ishing stress, (ii) the shear stress becomes negative at the 
end of the first cycle, (iii) no overshoot in the stress and 
the elastic energy from the second cycle, and (iv) that 
the stress and the elastic energy take roughly constant 
values characteristic of well-developed plastic flow in the 
region 0.25 < j(t) < 0.5 for increasing 7(2) and in the 
region < -f(t) < 0.25 for decreasing "f{t). Thus we can 
see significant hysteresis behavior. In MD simulations of 
low T glasses, similar stress-strain curves under step-wise 
strain rates have been obtained (but without overshoot 
behavior because of disordered initial states) [15,39]. 

We also notice that at the points A,B, and C, where 
the average stress vanishes as in Fig. 14, the curves of 
(/ e i) in Fig. 15 are locally minimum. This suggests that 
the strain j(t) consists of an elastic strain 7 e i and an slip 
strain 

7s = 7-7ei- (4.2) 

Roughly speaking, the elastic strain outside the disloca- 
tion cores should give to the average stress, while the slip 
strain is caused by the jumps of u x across the slips. To 
be more quantitative, we define 7 e i as 

4»(7ei) = (4-3) 

The elastic energy density stored is then the sum of the 
elastic energy density and the defect energy density. We 
define the average defect energy density by 

/b = </el>-*(7el,0), (4.4) 

where $ is defined by (2.5). For small 7 e i we have 
7ei — (<7 X y) and $(701,0) = 7ei/2. Hence, in the vicin- 
ity of the points where (a xy ) = 0, (f e \)(t) should take a 
minimum, provided that fr>(t) changes slowly there. To 
confirm these arguments, we plot fr>{t) in Fig. 16. For 
j(t) < 0.3 in the first cycle, however, it represents the 
elastic energy due to the inhomogeneous fluctuations of 
the local strains (mainly due to Se^) with the peak height 
at 0.012 (not shown in the figure). After this initial pe- 
riod, fr>{t) is in a range of 0.004 — 0.005 reasonably rep- 
resenting the elastic energy due to the defects produced 



in plastic flow. The energy variance J ((5f c \) 2 ) is also 
about 0.005 in plastic flow obviously due to the discrete 
nature of dislocation cores. For 7 = 10~ 4 , on the other 
hand, /b(i) is almost constant around 0.002 in plastic 
flow. 

B. Uniaxial stretching 

In Fig. 17 we show the normal stress vs strain at con- 
stant strain rate e for t > 0. The characteristic features 
are very similar to those in Fig. 10. At small e(t) = it 
these curves first follow the curve in the homogeneous 
case, 

2 

<J XX - <J yy = 2A 2 (e) = — [sin(Tre) + sin(27re)], (4.5) 

where A 2 was introduced in (2.15). For the curve of 
e = 10~ 4 and that of e = 10~ 3 (with the higher peak) 
we set u = at t = 0. They approach the elastic insta- 
bility point e = tt- 1 cos- 1 [(\/3^- l)/8] = 0.298. After 
the catastrophic formation of slips, the slip orientation 
angles are ±7r/4 with respect to the stretched direction 
in agreement with the experiments [7-10,27,34,35] and 
the simulations [16,17]. The discontinuity across the slip 
lines appears in the tetragonal strain e 2 , which can be 
seen in the two snapshots of Se 2 — e 2 — e. In Fig. 18 
the displacement vector u in a 1/4 region of the total 
system at e = 0.46 with e = 10~ 3 is displayed, where 
the orientations of the slips are the most favorable ones 
with ip — ±7r/4 as determined below (3.14). In Fig. 19 
four slips are placed at t = 0, exactly as in Fig. 13, and 
time-development of es, 5e 2 , and <5/ e i are followed. Here, 
however, these initial slip orientations are marginal in 
the sense that -F s ii P in (3.13) is independent of e. We find 
that a dislocation pair is newly created near each end of 
the preexisting slips to glide in the favorable directions 
with tp = ±7r/4. In the plastic flow regime at e = 0.4, 
they grow into thick layers containing high-density dislo- 
cations. 

Next we apply cyclic uniaxial stretching, where e(t) = 
10~ 3 in the time regions nt p <t< (n+l/2)t p and e(t) = 
— 10 -3 in the time regions (n + 1 /2)£ p < t < (n + l)t p . 
We choose t p = 1000, so the maximum of e(t) is 0.5 and 
the minimum is 0. Fig. 20 shows the stress-strain curve 
for the first four cycles, where it is nearly periodic from 
the second cycle. As in the shear deformation case we 
introduce the average elastic strain e e i by 

2A 2 (e cl ) = (a xx - a yy ) (4.6) 

and the average defect energy density /d by 

/D = (/ci)-$(0,e c i). (4.7) 

where $ is defined by (2.5). The average slip strain is 
given by e s = e — e e i. For small e e i we have 2e c \ = 



(pxx — &yy) and $(0, e e i) = Fig. 21 shows the av- 

erage elastic energy density (/ e i) and the defect energy 
density /d- The latter is in a range of 0.008 — 0.009 in 
plastic flow. The inset of Fig. 21 displays a snapshot of 
/ci at the point A in the first cycle, where the average 
normal stress difference vanishes. Comparing it with the 
snapshot of / e i in Fig. 15 under shear deformation, we no- 
tice a considerable difference in the spatial anisotropy of 
the dislocation distribution between the two cases. 

C. Strain-induced disordered states 

In the plastic flow regime dislocations are proliferated 
and a structurally disordered state is realized. This effect 
may be called strain-induced disordering. We mention 
a simulation by Ikeda et al. [55], who applied a tensile 
strain to a 3D model to induce a change from a perfect 
crystal to an amorphous solid. In Fig. 4 in our previous 
work [37] we switched off a shear flow (a) before the peak 
time of the stress, (b) just after the peak time, and (c) 
in well-developed plastic flow. Affine deformation in (a) 
was maintained, while no appreciable time evolution was 
detected after transients in (b) and (c). This means that 
the structurally disordered states are metastable. 

In Fig. 22, we show extremely slow relaxation of the 
average energy density (/ e i)(i), where a shear flow with 
7 = 10~ 3 is stopped at time £a — 600 (at the point A 
in Fig. 14) and the system is relaxed for a time period of 
t w = 10 5 . Here we set u x = jaLo and u y — at the top 
y = Lq with 7a = 0.40, together with u = at the bot- 
tom. The dislocation distributions here closely resemble 
that in the inset of Fig. 15 (are not identical because they 
are obtained from different runs). In the upper curve (a) 
the noise strength e t h in (2.25) is 0.1 as in the previous 
simulations in this section, where each maximum in the 
initial stage corresponds to an energy increase accom- 
panied with a configuration change around a dislocation 
core. From (a) the typical energy increase is of order 
10~ 5 A 2 <~ 0.1 for each event. However, there is no ap- 
preciable relaxation for t > 3 x 10 3 . In (b) we increase 
the noise strength e t h to 0.25 to obtain a larger energy 
decrease in (f e \)(t) with a larger thermal noise superim- 
posed. In these cases only a small number of configura- 
tion changes (~ 10) occur around dislocation cores even 
for t w — 10 5 , so the effect of the structural relaxation 
is negligible on the macroscopic level (no aging effect). 
In fact, the stress-strain curves after switching-on of the 
shear flow at t = + t w are almost independent of t w 
if plotted as a function of Aj(t) = A f(t — — t w ). In- 
terestingly, they exhibit a rounded peak at A7 = 0.15 as 
can be seen in the inset. In addition, in the waiting time 
region (£a < t < tA+t w ), the average (a xy )(t) fluctuates 
in time around 0, but its noise amplitude is only of order 
0.002 for (a) and 0.004 for (b) and the stress is nearly 
fixed. 



V. SUMMARY AND CONCLUDING REMARKS 



In summary, we have presented a nonlinear elasticity 
theory taking into account of periodicity of the elastic 
energy density with respect to the shear and tetragonal 
strains, e% and e-i. It has the symmetry of the 2D tri- 
angular lattice but is surprisingly isotropic in the e 3 -e2 
plane. We summarize our main results together with 
some comments. 

(i) We have numerically examined the slip structure as 
a function of its length £, its angle <p, and applied strains 
7 and e. In external strain field, slips should appear in 
the orientations minimizing the slip energy F s n p and they 
should grow into shear bands observed in previous ex- 
periments under external load in various materials. The 
snapshots of the displacement vector in the plastic flow 
regime, Fig. 12 for shear deformation and Fig. 18 for uni- 
axial stretching, most unambiguously illustrate the phys- 
ical processes taking place. We remark that the previous 
experiments have mostly been performed under uniaxial 
(biaxial or triaxial) deformation, so future experiments 
with shear deformation ({u x ) = jy and (u y ) = 0) should 
be informative. On the other hand, in crystalline solids 
with strong crystal anisotropy, the slip planes are paral- 
lel to particular crystal planes [56]. In our simulations 
we may obtain slips as steady solutions of our dynamic 
model due to the presence of the Peierls potential energy. 
However, they cannot be stationary for 7 < 7 c i <~ —0.1 
or 7 > 7 C 2 ~ 0.1 in shear strain 7, where the potential 
minima disappear and the slips shrink or grow. 

(ii) We have examined plastic flow by applying a con- 
stant strain rate at t = 0. If there is no initial disorder, 
the stress-strain curve exhibits a pronounced overshoot 
with a peak stress of order 10% of the shear modulus. 
However, in the presence of initial disorder, the over- 
shoot is weakened or even erased, as revealed by sim- 
ulations with initial four slips in Fig. 10, under cyclic 
shear in Figs. 14 and 20, and after staying at the zero- 
stress point A in the inset of Fig. 22. In accordance 
with these findings, previous simulations performed with 
various initial states have demonstrated sensitive depen- 
dence of the overshoot behavior on the quenching con- 
ditions [17,57]. In addition, a number of previous simu- 
lations have reported either of existence or nonexistence 
of the stress peak. In some real glassy systems including 
polymers, overshoot behavior has been widely observed 
[12,25-27]. In an amorphous metal, the stress increased 
monotonically to a steady-state value for slow strain rate 
(e < 10 _3 s _1 ), while a maximum appeared at e ~ 0.06 
for large strain rate (e > 5 x 10~ 3 s _1 ) [26]. 

(iii) As illustrated in Figs. 16 and 21 we have divided 
the total elastic energy into the affine part and the de- 
fect part following the definitions (4.4) and (4.7). In our 
model, once defects are created, the defect elastic energy 



is rather weakly dependent on the deformation history. 
On the basis of this division, we may easily understand 
the characteristic features of the stress-stain curves in 
the cyclic straining mentioned in the previous section. In 
future microscopic simulations on glassy materials, this 
kind of energy division should be informative at relatively 
small strains, where we are interested in the defect conri- 
bution to measurable quantities such as the specific heats. 
Here we mention microscopic calculations of the average 
potential energy per particle (e) in supercooled states un- 
der shear [57-59], where (e) was increased considerably 
by shear flow above the initial value in quiescent states. 
It is remarkable that, while (e) in quiescent states sen- 
sitively depends on the quenching history (aging effect) 
[57,60], it becomes uniquely determined in shear flow for 
shear rates larger than the inverse structural relaxation 
time r^ 1 . This is because sheared systems are effectively 
driven away from the glass transition [23]. This effect 
would be consistent with our result that /d(£) is kept 
nearly constant in plastic flow as in Figs. 16 and 21. 

(iv) In our simulations slips emerge as long straight 
lines, as shown in the snapshots of or e^. In our model 
the crystal order is not broken over long distances, but 
if disorder is fully introduced, the glide motions of slips 
in particular directions should be much limited. In MD 
simulations of two component glasses, for example, such 
degree of disorder should be sensitive to the size-ratio of 
the two species. In fact, it was rather close to 1 in the 
simulation by Deng et al. [12,15], where nano- crystalline 
order was realized [61] because of rather weak frustra- 
tion and long slips with atomic thickness emerged along 
the crystal axes. It is therefore informative to perform 
MD simulations with various size-ratios and examine the 
shapes of local configuration changes. 

(v) As a special ingredient of our theory, our elastic en- 
ergy is almost isotropic if the distance from the center of a 
unit cell in the es~e2 plane is shorter than 0.5 (see Figs.l- 
3). The snapshot at 7 = 0.4 in Fig. 11 demonstrates 
that most spatial points are in the mechanically stable 
regions in plastic flow and hence are in the isotropic elas- 
ticity regions. Therefore, our results such as the orienta- 
tions of well-developed shear bands should be applicable 
to those in amorphous materials (which are isotropic on 
large scales). 

(vi) The dimensionless parameter A in (2.23) is related 
to Poisson's ratio v as in (3.2) and the dislocation energy 
depends on v. Our choice Kq/hq = A = lor^ = 0is 
rather unusual, since K is usually considerably larger 
than /xo in high density systems. The appropriateness of 
the other choices t/q — 1 in (2.26) and e t h = 0.1 in Section 
4 should also be examined in future. 

(vii) As another aspect, we mention the effect of elas- 
tic interaction among dislocations. In our case slips are 
more easily created around preexisting ones, as already 
reported in Ref. [17]. In the insets of Figs. 15 and 21 we 
can see a tendency of aggregation of dislocation cores. 



In real 3D crystals a tangle of dislocations often appears 
as a characteristic feature of fatigued states [62] and has 
also been realized in 3D large-scale simulations [3,4]. In 
the literature the inter-dislocation elastic interaction is 
believed to yield mesoscopic patterns in the dislocation 
distribution [3,4,62]. 

(viii) As already stressed in the introduction, in order 
to describe glass dynamics in a more satisfactory level, 
we should try to include a variable representing the local 
free- volume [37] and the configurational frustration effect 
induced by the size disparity between the two species. 
This generalization should be essential with increasing 
T towards the glass transition. For example, we may 
predict a gradual diffusional increase of the free volume 
around dislocation cores [37], which induces configura- 
tion changes and presumably resulting in significant ag- 
ing effects. (The critical strains j c \ and j c2 discussed be- 
low (3.23) and indicated in Fig. 6 by arrows are sensitively 
decreased in magnitude by a small amount of the local 
free volume near the dislocation cores.) In the present 
study, as shown in Fig. 22, we have found no appreciable 
aging effect. 
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APPENDIX 

Here we examine the linearized equations of (2.16) 
and (2.17) around a homogeneously strained state with 
(ei) = 0, (e 3 ) = 7, and (e 2 ) = e. We neglect the ran- 
dom stress and assume that all the deviations depend on 
space and time as exp(ifc • r + iuot). Then the deviation 
of the displacement vector (Su x ,6u y ) obeys 

Q.5U X — C XX 5U X + C X ySUy, 

QSUy = C X ydU x + CyySUy, (A-l) 

where 

O = pu 2 /k 2 - ir) u>. (A. 2) 

The frequency u> is expressed in terms of at small k as 

w = ±(Q/p) 1/2 k + i{riv/2p)k 2 + 0(k 3 ). (A.3) 

The first term is the oscillation frequency for il > and 
the second term viscous damping. The coefficients C a p 
are expressed in terms of the coefficients $ Q( 3 in (2.7) as 



C xx = (K + $22)n 2 x + 2<& 23 n x n y + $33^, 

C yy = (K + <S>22)n 2 y - 2$23«x% + $33^, 

C xy = (Ko - $22 + * 33 K% + $23{n 2 x - n 2 y ). (A.4) 

Here n = k~ 1 k is the unit vector representing the di- 
rection of the wave vector k. From (A.l) it follows the 
relation, 

(n-c xx )(n-c yv ) = c 2 xy . (A.5) 

Let the angle of n be ip + tt/2 with respect to the x 
axis; then, 

n x = — sintp, n y — cosLp. (A. 6) 

After some calculations (A.5) is solved to give 

= 7^0 + $22 + $33) 

± K 2 + A 2 + B 2 + K ( A cos 4(^ + B sin Aip) , (A.7) 

where A = (<f> 22 - $33)/2 and B = $ 23 . As a func- 
tion of y the slowest mode is obtained if the combination 
Acos4< f 5+_Bsin4(^ takes the maximum {A 2 +B 2 ) 1 / 2 . The 
corresponding minimum of £1 is given by 

O min = \ ($22 + $33) - VA 2 +B 2 . (A.8) 

Notice that f2 m i n is the smaller of the two eigenvalues 
of the matrix {$ Q( g}. We may draw two conclusions, 
(i) As the instability point is approached, S! m i n tends to 
zero with softening of the sound speed (^min/p) 1 ^ 2 of the 
corresponding acoustic mode, (ii) If the approximate ex- 
pression (2.8) is used, we have A — 2G"(e 2 — -f 2 ) and 
B = 4G"e 7 so that 

Acos4tp + Bsin^ = Ccos(4<p + 2a), (A. 9) 

where C = 2\G"\(e 2 + 7 2 ) and a is defined by (3.15). 
Thus the minimum condition for f2 yielding (A.8) is given 
by 4(p + 2a = 2nir and is equivalent to (3.14). That is, 
the slowest mode, which undergoes softening at the in- 
stability point, has a wave vector perpendicular to the 
favorable slip orientations given by (3.14). This is the 
case even far below the instability point. In addition, for 
the slowest mode, (A.l) gives 

Su y /Su x = (f2 min - C xx )/C xy = cot ip, (A. 10) 

if a is eliminated using (3.14). Thus the deviation 5u is 
perpendicular to the wavevector k or the slowest mode 
is a transverse sound. 
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FIG. 1. The scaling function 3>(e3, ei) in (2.5), which is the shear deformation energy density divided by fio with one of the 
crystal axes being along the x axis. 



FIG. 2. $(e3,e2) for various directions in the ez — ei plane. It demonstrates isotropic behavior for e = (e| + e%) 1 ^ 2 < 0.5. 
The curve for uniaxial stretching (ei = e and e$ — 0) coincides with the curve of n = 5 in the figure. The instability points are 
marked, which separate stable and unstable regions 



FIG. 3. The stable regions of our elastic energy are shown in white, where the two eigenvalues of {$ Q /3} are both positive. 
In the gray region one of them is negative, and in the black regions both of them are negative. 



FIG. 4. The x component of the displacement for a numerically calculated type C slip with length 20 along the x axis in the 
unstrained condition. 



FIG. 5. The displacement vector u around type C and type CC slips with length 10. The orientations in (a) are most 
favorable in shear deformation, while those in (b) are most favorable in uniaxial stretching. The arrows are from the original 
undeformed position to the displaced position. 



FIG. 6. The slip energy F s n p vs applied shear strain 7 for type C slips with length i — 10, 20, and 30 oriented along the 
x axis. For | -7- 1 < 0.05 the relation (3.12) holds with slope —I. The arrows indicate the instability points where expansion or 
shrinking of the slips occur. 



FIG. 7. The slip energy density around a type C slip with length 20 for 7 = in (a), 0.065 in (b), and -0.04 in (c). In the 
middle region between the dislocations at the ends, the elastic energy is decreased in (b) and increased in (c) . 



FIG. 8. Coincidence of the derivative dF B \ ip /d~/ (x) and the space integral of a xy — 7 (+) confirming the general relation 
(3.17). These quantities are calculated from the steady slip solutions. 



FIG. 9. 

FIG.9. The slip energy difference AF slip {£) = F slip (£) - F s r lp (20) in the range 20 < I < 21 obtained for the 
extrapolated displacement (3.22) for 7 = —0.04,0,0.08, and 0.1. The position of the maximum is a decreasing 
function of 7. The maximum position approaches £ = 21 as 7 — > 7 c i = —0.09 and £ — 20 as 7 — ► j C 2 = 0.12. The 
inset shows the displacement vector it 2 i — 1*20 needed for slip growth by unit length from £ = 20 to 21. 

FIG. 10. The stress-strain curves in shear flow. The initial states are defectless for 7 = 10~ 4 and 10~ 3 (with the higher 
peak). Four slips are initially prepared for 7 = 10~ 3 (with the lower peak). Embedded are snapshots of <5e3 = e-z — 7 at 7 = 0.3 
at the inception of slip formation and at 7 = 0.4 in the plastic flow regime. 
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FIG. 11. Mechanically stable regions (white) and unstable points (black) after application of shear flow with 7 = 1CP 3 . The 
initial state is crystalline without defects. At the stable points the eigenvalues of the matix {$ a /3} defined below (2.7) are both 
positive, while at the unstable points at least one of them is negative. At 7 = 0.29 the system is close to the peak position, at 
7 = 0.30 slips are appearing, and at 7 = 0.4 the unstable points are fluctuating in time near the dislocation cores. 



FIG. 12. The displacement deviation 5u = u — (73/, 0) in the plastic flow regime under shear strain at 7 = 0.4 with 7 = 10~ 4 . 
A 1/4 region (64 x 64) of the total system is shown. The arrows are from the original position at t — in a perfect crystal to 
the displaced position in plastic flow. 



FIG. 13. Time-evolution of <5e3 and 5f e \ in shear flow when four slips are placed in the initial state. 



FIG. 14. 

FIG. 14. The stress-strain curve for cyclic shear deformation in the first two cycles at 7 = ±10~ 3 with period 10 3 . 
Once 7 > 0.3 in the first cycle, high-density dislocations are created. The shear stress vanishes at the three points A, 
B, and C (x). 

FIG. 15. The average elastic energy density (/ e i) vs the strain 7(f) for cyclic shear deformation in the first two cycles. At 
the points A, B, and C the shear stress vanishes. In the inset the snapshot of f e \ at the point A is shown, where (a xy ) = 7 e i = 
and the black points represent dislocation cores. 



FIG. 16. The energy density /d defined by (4.2) vs the strain -y(t) for cyclic shear deformation in the first two cycles. It 
is the defect energy density in plastic flow, but in the preplastic regime it arises from the heterogeneities in the strains and is 
enhanced at the onset of plastic flow. 



FIG. 17. The stress-strain curves under uniaxial stretching. The initial states are defectless for e = 10~ 4 and 10~ 3 (with the 
higher peak). Four slips are initially prepared for e = 10~ 3 (with the lower peak). Embedded are snapshots of Se? — e2 — e at 
e = 0.31 at the inception of slip formation and at e — 0.46 in the plastic flow regime. 



FIG. 18. The displacement deviation Su = u — (ex/2, —ey/2) in the plastic flow regime under uniaxial stretching at e = 0.46 
with £ = 10~ 3 . A 1/4 region (64 x 64) of the total system is shown. The arrows are from the original position at t — in a 
perfect crystal to the displaced position in plastic flow. 



FIG. 19. Time-evolution of es, 5e2, and 8f e \ under uniaxial stretching when four slips are placed at t = 0. The initial 
orientations are marginal and new slips with the favorable orientations appear from the ends of the initial slips. 



FIG. 20. 

FIG. 20. The stress-strain curve for cyclic uniaxial stretching in the first four cycles at e = ±10~ 3 with period 10 3 . 
Once e > 0.3 in the first cycle, high-density dislocations are created. In this case the stress-strain relation becomes 
nearly periodic from the second cycle. 
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FIG. 21. The average total elastic energy density (/ e i) (solid line) and the defect part fo (broken line) vs the strain e(t) for 
cyclic uniaxial stretching in the first two cycles. In the inset the snapshot of / c i at the point A in the first cycle in Fig. 20 is 
shown, where (a xx — o yy ) — e c i = and the black points represent dislocation cores. 



FIG. 22. Relaxation of the average energy density (/ e i)(i) for eth = 0.1 (a) (upper curve) and eth = 0.25 (b) (lower curve). 
The shear flow is stopped at time £a = 600 (at the point A in Fig. 14) and the system is relaxed for a time period of t w = 10 5 . 
The energy decreases are extremely small as compared to the initial values, so there are almost no appreciable aging behavior 
here. The inset displays the stress-strain curves, (a X y) vs A7(t), after the shear is switched on again for t > £a + iw, where the 
solid line corresponds to (a) and the noisy broken line to (b). The A7(t) = j(t — £a — t w ) is the excess strain with 7 = 10~ 3 . 
Almost identical stress-strain curves are obtained for any waiting time t w shorter than 10 5 . 
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